home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
cpp_libs
/
newmat03.lha
/
newmat03
/
hholder.cxx
< prev
next >
Wrap
C/C++ Source or Header
|
1993-08-08
|
1KB
|
56 lines
//$$ hholder.cxx Householder triangularisation
// Copyright (C) 1991: R B Davies and DSIR
#define WANT_MATH
#include "include.hxx"
#include "newmatap.hxx"
/********************** householder triangularisation *********************/
inline real square(real x) { return x*x; }
void HHDecompose(Matrix& X, LowerTriangularMatrix& L)
{
int n = X.Ncols(); int s = X.Nrows(); L.ReDimension(s);
real* xi = X.Store(); int k;
for (int i=0; i<s; i++)
{
real sum = 0.0;
real* xi0=xi; k=n; while(k--) { sum += square(*xi++); }
sum = sqrt(sum);
L.element(i,i) = sum;
real* xj0=xi0; k=n; while(k--) { *xj0++ /= sum; }
for (int j=i+1; j<s; j++)
{
sum=0.0;
xi=xi0; real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
L.element(j,i) = sum;
}
}
}
void HHDecompose(const Matrix& X, Matrix& Y, Matrix& M)
{
int n = X.Ncols(); int s = X.Nrows(); int t = Y.Nrows();
if (Y.Ncols() != n) MatrixError("Incompatible dimensions in HHDecompose");
M.ReDimension(t,s);
real* xi = X.Store(); int k;
for (int i=0; i<s; i++)
{
real* xj0 = Y.Store(); real* xi0 = xi;
for (int j=0; j<t; j++)
{
real sum=0.0;
xi=xi0; real* xj=xj0; k=n; while(k--) { sum += *xi++ * *xj++; }
xi=xi0; k=n; while(k--) { *xj0++ -= sum * *xi++; }
M.element(j,i) = sum;
}
}
}